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Summary 

The Integrated Force Method, a recently developed method 
for analyzing structures, is extended in this paper to three- 
dimensional structural analysis. First, a general formulation 
is developed to generate the stress interpolation matrix in terms 
of complete polynomials of the required order. The formula- 
tion is based on definitions of the stress tensor components in 
terms of stress functions. The stress functions are written as 
complete polynomials and substituted into expressions for 
stress components. Then elimination of the dependent coeffi- 
cients leaves the stress components expressed as complete 
polynomials whose coefficients are defined as generalized 
independent forces. Such derived components of the stress 
tensor identically satisfy homogenous Navier equations of 
equilibrium. The resulting element matrices are invariant with 
respect to coordinate transformation and are free of spurious 
zero-energy modes. The formulation provides a rational way 
to calculate the exact number of independent forces neces- 
sary to arrive at an approximation of the required order for 
complete polynomials. The influence of reducing the num- 
ber of independent forces on the accuracy of the response is 
also analyzed. The stress fields derived are used to develop a 
comprehensive finite element library for three-dimensional 
structural analysis by the Integrated Force Method. Both 
tetrahedral- and hexahedral-shaped elements capable of mod- 
eling arbitrary geometric configurations are developed. A 
number of examples with known analytical solutions are 
solved by using the developments presented herein. The re- 
sults are in good agreement with the analytical solutions. The 
responses obtained with the Integrated Force Method are also 
compared with those generated by the standard displacement 
method. In most cases, the performance of the Integrated Force 
Method is better overall. 


Introduction 

The finite element method has become the preferred tool 
for analyzing a wide variety of physical problems. The 
advantages of the finite element method over other numerical 
techniques, such as the finite difference method and the bound- 
ary element method, include the efficient and accurate mod- 
eling of domains with arbitrary geometric configurations and 
varying material parameters and the capability to analyze prob- 
lems with both material and geometric nonlinearities. Two 
distinct finite element formulations, based on the primary 
unknown variables used in the analysis, have been developed 
for analyzing problems in structural mechanics, namely, the 
displacement method (refs. 1 to 3) and the standard force 
method (refs. 4 and 5). The displacement method has become 
prevalent in structural analysis because of its straightforward 
implementation and efficient use of computer resources. It 
has been implemented in all commercial finite element pro- 
grams. Implementation of the standard force method, on the 
other hand, requires either a selection of redundant forces 
(ref. 4) or extensive manipulation of system matrices (refs. 5 
to 7) in order to generate a system of equations for the un- 
known variables. Since these procedures could not easily be 
adapted for computer automation and since they lacked a 
physical interpretation, the standard force method met its 
demise. 

Certain drawbacks of the displacement-based method have 
been observed in a number of applications, such as the bend- 
ing of thin plates, the analysis of nearly incompressible mate- 
rials, and the optimization of structures (refs. 7 to 11). In the 
displacement method, stresses are calculated indirectly by 
using displacement derivatives, which may introduce errors 
in stress predictions. Thus, there is a need to develop alterna- 
tive finite element formulations to treat the aforementioned 
problems and to provide a means of comparison for other 
numerical methods. 



A new force method formulation for analyzing structures 
(refs. 12 to 15), known as the Integrated Force Method, has 
been developed in recent years. In the Integrated Force Method 
all independent forces are treated as unknown variables that 
can be obtained by solving a system of equations consisting 
of equations of equilibrium and compatibility conditions. 
Compatibility conditions are derived in a procedure similar 
to the St. Venant procedure for continuous elasticity (ref. 16). 
Several classes of compatibility conditions have been identi- 
fied, and a method for generating compatibility conditions 
has been developed (ref. 17). This method efficiently utilizes 
computer resources and produces a sparse and banded com- 
patibility matrix. The Integrated Force Method has been suc- 
cessfully applied to static (ref. 12) and free vibration (ref. 18) 
analyses and to structural optimization (refs. 19 to 21). Initial 
comparisons with other finite element formulations have re- 
vealed the superiority of the Integrated Force Method in ac- 
curacy as well as in computer efficiency (ref. 22). 

In the Integrated Force Method, two sets of relations are 
established: the equilibrium relations and the deformation- 
force relations. These relations are first derived on the ele- 
ment level through equilibrium and flexibility matrices, 
respectively (ref. 23). Next, appropriate assembly techniques 
are used (ref. 14) to generate systems of equations for the 
entire structure. In order to generate element matrices, approxi- 
mations of both displacement and stress fields have to be 
defined. Here, these approximations are performed indepen- 
dently; the displacement components within an element are 
approximated in terms of element nodal displacements, 
whereas the components of the stress tensor are proposed in 
terms of a set of independent generalized forces. Although 
the interpolation of the displacement field is straightforward 
and does not significantly influence the accuracy of calcula- 
tions in the Integrated Force Method, the correct approxi- 
mation of the stress field is of utmost importance for both 
accuracy and computer efficiency. Stress field interpolations 
were previously studied in applications of the hybrid method 
(refs. 8, 9, and 24 to 26). Spilker and Singh (ref. 9) derived 
two stress field interpolations for a 20-node isoparametric 
element; Pian and Sumihara (ref. 25) used nondimensional 
coordinates; and Pian et al. (ref. 26) explored the possibility 
of achieving the minimum number of independent forces 
(ref. 27) while relaxing some requirements for stress fields. 

The Integrated Force Method is extended in this paper to 
three-dimensional structural analysis. To this end an exten- 
sive library of finite elements capable of analyzing domains 
with arbitrary geometric configurations is developed. First, a 
general formulation is developed to generate the stress inter- 
polation functions in terms of complete polynomials of the 
required order. The formulation is based on definitions of stress 
components in terms of stress functions. The stress functions 
are written in terms of complete polynomials of a certain or- 
der with coefficients that are arbitrary constants. Expressions 


for stresses are obtained next by substituting expressions for 
stress functions. Then the linearly dependent constants are 
eliminated, thereby yielding final stresses expressed in terms 
of complete polynomials whose coefficients are now element- 
generalized independent forces. Such derived stress fields 
identically satisfy Navier’s equations of equilibrium. The el- 
ement matrices generated with these polynomials are invari- 
ant with respect to coordinate transformation and are free of 
spurious zero-energy modes for arbitrary orientation of 
the element local coordinate system. This formulation also 
provides a rational way to calculate the exact number of inde- 
pendent forces for an approximation with complete polyno- 
mials of the required order. Representing the stress field with 
complete polynomials may result in a significantly larger 
number of independent forces than the minimum number 
determined by rigid body modes considerations. This paper 
discusses the reduction in the number of independent forces 
and some guidelines for determining the stress fields proper- 
ties necessary to achieving good approximations. 

The stress fields derived by using the developments pre- 
sented herein are used to develop a comprehensive finite ele- 
ment library for three-dimensional structural analysis by the 
Integrated Force Method. Both tetrahedral- and hexahedral- 
shaped elements capable of modeling arbitrary configurations 
are' considered. Stress fields that were derived earlier for 
hybrid elements (refs. 5 and 9) and given in terms of reduced 
polynomials are also implemented to study the effects of 
reduction on accuracy. 

To establish their validity and accuracy, the elements pre- 
sented here are used to solve a number of example problems 
with known analytical solutions. These results are compared 
to those obtained with the standard displacement method to 
assess the relative performances of the two methods. The ex- 
pressions for stress fields are given in appendix A, and the 
symbols used are defined in appendix B. 

Basic Equations of the Integrated 
Force Method 

The governing equations of the Integrated Force Method 
are briefly reviewed here to ensure completeness and to in- 
troduce the notation. The derivation of these equations and a 
description of the method are detailed in references 12 to 14. 
Expressions for element matrices are derived in reference 23. 

In an analysis by the Integrated Force Method, a continu- 
ous object under consideration is discretized into a number of 
finite elements. Within each finite element the displacement 
and stress fields are approximated in terms of two sets of in- 
dependent variables as 

{u} = [N]{U e } (1) 
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{t7} = [Y]{FJ (2) 

where {u} T = {u v w } is the displacement vector with com- 
ponents w, v, and w along the global coordinate axes Ox , Oy , 
and Oz , respectively; {U e } is the vector of displacements at 
element nodes resulting from the finite element discretization; 
{a} T = {<y x a y g z X xy x yz x^ is the vector of stress compo- 
nents; {F e } is the vector of element independent generalized 
forces; [N] is the matrix of displacement interpolation func- 
tions; and [Y] is the stress interpolation matrix. Such a finite 
element discretization results in a total of N t displacement 
degrees of freedom and m independent forces. Two sets of 
relations can be written for each element of the discretization: 

(a) the equations of equilibrium 

{P c } = [BJ{F e } (3) 


and 


be augmented by a set of r = m - n compatibility conditions 
(ref. 17) which, for the case of mechanical loads, have the form 

[C][G]{F} = {0} (7) 

where [G] is the system flexibility matrix and [C] is the 
compatibility matrix. The procedures for automatically gen- 
erating the compatibility matrix are given in reference 17. 
Equations (6) and (7) provide the necessary and sufficient 
number of equations to calculate the vector of independent 
forces {F}. The resulting system of equations may be written 
in a compact form as 

[S] {F} = {P*} (8) 

where 


[S] = 


' CB p ] ' 
JC][G]_ 


and {P*} = | { { J } } } 


( 9 ) 


(b) the deformation-force relation 

{P € } = [GJ{F € } (4) 

where {P e } is the equivalent nodal load vector of the element; 
{P g } is the vector of element deformations corresponding to 
the forces {F^}; and [B e ] and [G e ] are the element equilib- 
rium and flexibility matrices, respectively, which are calcu- 
lated (ref. 23) as 

[B e ] = f [Z] T [Y]dV (5a) 

Jv 


After the force vector is calculated from equation (8), the vec- 
tors of the unknown nodal displacements {U} and the sup- 
port reactions {R} may be obtained as 

{U} = [J] [G] {F} (10a) 

{R} = [BJ {F} (10b) 

where [ J] is the n x m deformation matrix that represents the 
top m rows of the transpose of the matrix [S]*** 1 ; and [ B s ] is 
the portion of the system equilibrium matrix that corresponds 
to nodes with prescribed displacement boundary conditions. 


[GJ = J v [Y] T [D][Y]dV (5b) 

In equations (5), [D] is the compliance matrix of the material, 
and [Z] = [L][N], where [L] is the matrix of differential 
operators that defines the strain-displacement relationship. For 
all free nodes, equation (3) can be written to generate a sys- 
tem of n equations of equilibrium: 

[Bp]{F} = {P} (6) 

Here [B p ] is the {m x n) system equilibrium matrix; {F} is 
the system vector of independent forces of dimension n; and 
{ P } is the m-component vector of equivalent nodal loads. For 
a general problem, m > n, and the system given in equa- 
tion (6) is indeterminate and not sufficient to calculate the 
independent forces. This system of equations (eq. (6)) must 


Approximations of the Stress Field 

Approximating the stress field within a finite element and 
generating the stress interpolation matrix [Y] are discussed in 
this section. The stress interpolation matrix appears in defini- 
tions of both the equilibrium and flexibility matrices (see 
eqs. (5)). It is, therefore, necessary to properly devise the stress 
interpolation functions in order to obtain an accurate response. 
The required properties of stress interpolation polynomials, 
in the context of hybrid method applications, were discussed 
by Spilker et al. (refs. 8 and 9) and Pian et al. (ref. 26). They 
found that approximate stress fields should identically satisfy 
Navier’s equations of equilibrium, that the stress components 
should be symmetric, and that the resulting element matrices 
should be invariant with respect to coordinate transformations 
and be free of spurious zero-energy modes. Spilker, Masked, 
and Kania (ref. 8) have shown that a necessary and sufficient 


3 



condition for element matrices to be invariant with respect to 
coordinate transformation is that stress fields be approximated 
in terms of complete polynomials. In reference 23 , the same 
requirement was deduced for the Integrated Force Method fi- 
nite elements. A formulation to derive stress interpolation 
functions as complete polynomials of arbitrary order is de- 
veloped herein, and then the reduced stress fields are dis- 
cussed. 

Stress Fields in Terms of Complete Polynomials 

Stress interpolation functions in terms of complete poly- 
nomials of order p- 2, where p > 2 is an integer, are derived 
in this section. To this end, the stress functions fork- 1 , 
2 , 3 , are proposed as complete polynomials of order p such 
that 

p P~i 

® k (x,y,z) = '^^(f)x i y i z p - i - j k = 1,2 , 3 (11) 
1=0 7=0 


Shape 

Nodes 

Forces 

Name 

Number of 
integration 
points 


4 

6 

THD04JJ6 

1 

18 

THD04J8 

5 

21 

THD04_21 

5 

: 1£ 

8 

18 

HEX08.18 

3x3x3 

* 


33 

HEX08_33 

3x3x3 

48 

HEX08_48 

3x3x3 

VA 


10 

36 

THD1 0_36 

11 

39 

THD1 0_39 

11 

48 

THD10_48 

11 



20 

57 

HEX20_57 

4x4x4 

60 . 

HEX20__60 

4x4x4 


90 

HEX20_90 | 

4x4x4 


Figure 1 . — Finite elements for three-dimensional analysis. 
(Numerical integration for THD elements use ref. 1; likewise, 
HEX elements use ref. 29 for quadrature.) 


where Cfj for i - 0, 1, ..., p, and j = 0, 1, ..., ip - i) are arbi- 
trary coefficients, and (x, y, z) denotes the position of a point 
within a finite element in the element local coordinate sys- 
tem. Local coordinate systems for various element shapes are 
depicted in figure 1 . The expressions for stress components 
can be obtained from the definitions of stress functions 
(ref. 28 ) as 




i =0 7=0 

+ + 1)() + 2>] ^-W-2 (12a) 


p-2 p-i —2 


Cy dx 2 


i =0 7=0 

+ Cg> (p - i -j){p - i - j - 1)] x i y j z p ~ H ~ 2 (12b) 


p —2 p-i —2 


i =0 7=0 

+ C[ 2 lj(i + 1 )(i + 2)] x i y j z p ~ i ' j ~ 2 (12c) 


_£_L. 

^ dxdy 


i 2 >k p zlz 2 


— X X 


i= 0 7=0 


x x i y i z p ~ i ~ i ~ 2 (12d) 


i ~0 7=0 


x x‘y J z p 1 7 2 (12e) 


_ d 

%zx " dzdx 


p -2 p-i -2 




i=0 j =0 


x X i y j z p ~ i ~ j ~ 2 ( 12 f) 


Equations ( 12 ) can now be rewritten as 


p—2 p-i —2 


Iv 


i =0 7=0 


p-2 p-i —2 

-XX 


^ X i y i z p ~ i ~ j ~ 2 
■2 )+/ * y z 




( 13 a) 


i=0 y=0 


( 13 b) 
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p-2 p-i- 2 

G ; = X X f 2P + ix{p-i-l) + i ^ i y^ p ~ i ~ hl (13c) 

(=0 7=0 

p- 2 p-i- 2 

s=X ° 3,i) 

i=0 7=0 

p-2 p-i- 2 

x yz = X X P 4P + ix(p-i-,) +j x i y 3 Z p ~ H ~ 2 (13e) 

r =0 7=0 

p-2 p-i- 2 

v - X X WiM-w <13f) 

1=0 7=0 

where P = (l/2)/?(p - 1), and the coefficients F ? , for # = 1, 
2, 6P, denote element-generalized forces. Generalized 

forces F q are expressed in terms of C® as 

f ff M*(cg) for q = 1,2, ...,6P (14) 

where <t>^ represents the linear functions of constants C®. 
From equation (14) we can see that 3p(p - 1) forces F q are 
expressed in terms of (3/2 )(p - l)(p + 4) constants C®. 
Eliminating all dependent quantities from equations (13) yields 
final expressions for stress polynomials in terms of Q = (3/2) 
x (p - l)(p + 2) independent generalized forces F qi for q = 1, 
2, 2. Such generated stress fields identically satisfy the 

equations of equilibrium. The element matrices generated by 
using these stress fields are invariant with respect to coordi- 
nate transformation and are free of spurious zero-energy 
modes. 

A general procedure for deriving the stress interpolation 
matrix in terms of complete polynomials can be illustrated 
through the derivation of linear terms. For this case, the stress 
functions <£>*., for k = 1, 2, 3, are written as complete third- 
order polynomials. 


components in terms of stress functions, the expressions for 


stresses can be obtained as 


= («g + 2Cg)*+(2Cg + 2Cgl)y 


+(6Cgo> + 2Cg) z 

(16a) 

= ( 6C S + 2C$)* + (2Cg + 2Cg )y 


+(2Cg>+6Cg>) Z 

(16b) 

0 Z = (2Cg + 2Cg>)x + (6Cg> + 2CfJ)y 


+ (2C$ + 2Cg>) Z 

(16c) 


(16d) 

vH3V* +2 ^ + 2Cg> Z ) 

(16e) 

V = -( 2C S^ + C U>' + 2Cg Z ) 

(16f) 

Equations (16) may be written in terms of element forces F t 
as 

g x = FjX + F 2 y + F 3 z (17a) 

o y = F 4 x + F 5 y + F 6 z 

(17b) 

a 2 =F 7 x+F i y+F 9 z 

(17c) 

x xy = ho x+ hiy + ht z 

(17d) 


® k {x,y,z) = C$ z 3 + C$yz 2 + C$y 2 z+C$y 3 + C$z 2 x 


= F u x+F u y+F l5 z 


+ C[ k ^xyz + C[ k ^xy 2 + x 2 z + C^x 2 y + cgjjjt 3 (15) 

where d k> , for i = 0, 1, 2, 3, and j = 0, 1, .... (3 - i) are coeffi- 
cients of tiie cubic polynomials. From the definitions of stress 


x z X = Fi6 x+ hyy + hs z ( 17f ) 

Expressions for forces F t in terms of coefficients can 
be obtained by comparing corresponding terms in equations 
(16) and (17). A careful examination of these expressions 
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reveals that not all forces F/ are linearly independent. Three 
relationships exist: 

Fi + Fn + Fig = 0 (18a) 

^ + ^10 + ^15 =0 (18b) 

(18°) 

Eliminating forces Ffy, F 15 , and Fjg in equations (17) 
through use of equations ( 1 8) and renumbering yields the lin- 
ear terms of the stress polynomials given in equations (A2) in 
appendix A. The terms of constant, quadratic, and cubic or- 
ders can be obtained by following a similar procedure; they 
are given in equations (A 1), (A3), and (A4), respectively. The 
stress field representation in terms of complete polynomials 
of order p can be obtained by combining the expressions for 
orders 0, 1 The interpolation in terms of complete cubic 
polynomials is given in equations (A5). The constant stress 
field is obtained by retaining the first six terms in equa- 
tions (A 5). For the linear interpolation, the first 21 terms are 
retained; for the quadratic interpolation, 48 terms are retained. 

Reduced Stress Fields 

Stress interpolation with complete polynomials may result 
in a large number of independent forces for each element, 
which leads to a final system of equations that is quite large. 
This problem is particularly pronounced in three-dimensional 
analyses, where the difference between the number of inde- 
pendent forces in complete polynomials and the minimum 
required number calculated from the number of rigid body 
modes for a particular element (ref. 27) grows rapidly as the 
order of interpolation increases. An effort should therefore be 
made to reduce the number of independent forces in stress 
interpolation polynomials while preserving the accuracy and 
reliability of the resulting elements. Reduced stress fields were 
studied by a number of researchers (refs. 9 and 26) with the 
goal of devising stress interpolations containing the minimum 
number of independent forces. Such derived stress fields, how- 
ever, may violate the requirements for accuracy stated earlier. 
There is no rational procedure available to uniquely 
derive stress fields with the minimum number of independent 
forces, and there is no proof that the resulting elements are 
free of spurious zero-energy modes. Moreover, in some prob- 
lems these elements may fail unexpectedly. Therefore, reduced 
stress fields should be used with caution, and extensive nu- 
merical studies should be performed to verify the resulting 
elements. The guidelines suggested by Pian, Chen, and Kang 
(ref. 26) were followed in this study to derive the two reduced 
quadratic stress fields given in equations (A9) and (AlO) and 
the reduced cubic stress field given in equations (All). The 


stress fields developed by Spilker and Singh (ref. 9) and 
Robinson (ref. 5) for hybrid finite elements were also imple- 
mented for comparison purposes. 

Finite Elements for Three-Dimensional 
Analysis 

In this section, a comprehensive finite element library for 
three-dimensional analysis by the Integrated Force Method is 
presented. Both tetrahedral- and hexahedral-shaped elements 
capable of modeling domains with arbitrary configurations 
were generated. Two groups of elements were developed for 
each shape. In the first group, nodes were introduced only at 
the vertices, thereby producing four-node tetrahedrons and 
eight-node hexahedrons. In the second group, additional nodes 
were introduced on the element edges as well, thereby pro- 
ducing 10-node tetrahedrons and 20-node hexahedrons. The 
attributes of the elements of the library presented here are 
depicted in figure 1. The names given to these elements con- 
sist of three parts: the first three characters denote the ele- 
ment shape, two subsequent digits denote the number of 
element nodes, and the number following the underscore de- 
notes the number of independent forces used in the stress 
interpolation. A local coordinate system Oxyz for stress inter- 
polation is defined such that the origin O coincides with the 
element centroid, and the local axes Ox, Oy , and Oz are paral- 
lel to the corresponding global axes. For all elements presented 
here, isoparametric functions (ref. 1 ) are used to inteipolate 
the displacement fields. The characteristics of these elements 
are enumerated in the following sections. 

Four-Node Tetrahedral Elements: THDO4_06, 
THD04JL8, and THD04_21 

Four-node tetrahedrons have 12 kinematic degrees of free- 
dom and thus require at least 6 independent forces in their 
stress field description. The six independent forces define com- 
plete polynomials of order zero, that is, the constant stress 
field. This stress field was implemented for element 
THD04J36. Since element THD04_06 contains the minimum 
required number of independent forces, it is a statically deter- 
minate element (ref. 23). Higher order stress fields were also 
implemented to investigate the influence that the order of in- 
terpolation has on accuracy. The complete linear stress field 
was used for element THD04_21 , and the 18-force stress field 
developed by Robinson (ref. 5), given in equations (A 6 ), was 
used for element THD04J 8 . 

Eight-Node Hexahedral Elements: HEX08_18, 
HEX08.33, and HEX08_48 

Eight-node hexahedral elements have 24 kinematic degrees 
of freedom and thus require at least 18 independent forces in 
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their stress field interpolation polynomials. Robinson’s stress 
field, given in equations (A6), was implemented for element 
HEX08_18, which represents a statically determinate eight- 
node hexahedron. Let us consider the complete polynomials 
for eight-node hexahedrons. First-order polynomials with 21 
independent forces satisfy the criterion for the minimum num- 
ber of forces. The test for zero-energy modes (ref. 23) reveals, 
however, that the resulting element possesses three spurious 
zero-energy modes. This behavior confirms that displacement 
and stress fields within an element cannot be chosen arbitrarily, 
but should be compatible with the stress-strain relationships. 
In order to obtain an element with complete polynomials, die 
second order was employed for element HEX08_48. Since 
this element contains a significantly larger number of inde- 
pendent forces than necessary, a reduced stress field was ob- 
tained with 21 independent forces corresponding to complete 
linear polynomials, combined with a reduced quadratic field 
with 12 additional forces (given in eqs. (A8)). This stress field 
identically satisfies Navier’s equations of equilibrium, and the 
resulting element does not possess spurious zero-energy 
modes. 

Ten-Node Tetrahedral Elements: THD10_36, THD10„39, 
and THD10_48 

Ten-node elements possess 30 kinematic degrees of free- 
dom and require 24 independent forces in stress interpolation 
polynomials. For this case, at least a second-order polyno- 
mial is required, and such an implementation was carried 
out for element THD10_48. Two reduced quadratic stress 
fields were also implemented. Both of these fields contain 
complete linear polynomials and a number of quadratic terms. 
The linear polynomials were combined with the quadratic 
terms given in equations (A9) to generate the stress fields used 
for element THD10J39, and with the quadratic terms given 
in equations (A10) to obtain the stress fields implemented for 
element THD10J36. Both of these fields identically satisfy 
the equations of equilibrium, and the resulting elements are 
free of spurious zero-energy modes. 

Twenty-Node Hexahedral Elements: HEX2CL57, 
HEX20J50, and HEX20_90 

Twenty-node hexahedral elements have 60 kinematic de- 
grees of freedom; thus, at least 54 independent forces must 
be present in their stress interpolation polynomials. Second- 
order complete polynomials contain 48 independent forces — 
not enough for the 20-node elements— so cubic polynomials 
must be used. Complete third-order polynomials were imple- 
mented for element HEX20_90, and the stress field devel- 
oped by Spilker and Singh (ref. 9) was used for element 
HEX20_57. An additional reduced stress field with 60 inde- 
pendent forces was developed in this study and implemented 


for element HEX20_60. It contains 48 independent forces 
corresponding to complete quadratic polynomials, and 12 ad- 
ditional forces representing cubic terms (given in eqs. (All)). 
Both reduced stress fields identically satisfy the equations of 
equilibrium, and neither element HEX20_57 nor element 
HEX20JS0 possesses zero-energy modes. 

For all elements developed here, numerical integration was 
used to calculate element matrices. For the tetrahedral ele- 
ments, a 1 -point rule that employs the constant stress field 
was used for element THD04_06; a 5-point rule was used for 
element THD04_21; and an 11-point rule was used for all other 
tetrahedral elements. The locations for integration points and 
the corresponding weights were taken from reference 29. Stan- 
dard Gauss integration was used for hexahedral elements, with 
3x3x3 points for elements with quadratic stress interpola- 
tions, and 4x4x4 points for elements with cubic stress 
interpolations. 


Example Problems 

A number of example problems are presented in this sec- 
tion in order to demonstrate the accuracy and validity of ele- 
ments developed in this study. Problems from beam theory, 
plane stress, and plate bending, for which analytical solutions 
are available, were selected. Extensive numerical experiments 
were performed to assess relative performances of the present 
elements and to compare the Integrated Force Method with 
the standard displacement method. The responses for one- 
dimensional problems analyzed here with three-dimensional 
discretizations are compared to those obtained with two- 
dimensional models (ref. 23), to verify analogous behavior in 
corresponding elements. 

Example 1: Bending of a Uniform Cantilever Beam 

A cantilever beam of length L with a uniform rectangular 
cross section of dimensions d by H is shown in figure 2(a). 
The beam is assumed to be made of a homogeneous and iso- 
tropic material with a modulus of elasticity E and a Poisson’s 
ratio v; it is subjected to a concentrated force of intensity P at 
the free end. This beam was analyzed with the entire finite 
element library in order to verify the present elements and to 
assess their relative performances. The influence of element 
shapes on the results was also analyzed. A three-dimensional 
finite element discretization of the beam is shown in fig- 
ure 2(b). The support conditions at the clamped end were 
modeled so as to suppress all three components of the dis- 
placement at point a as well as to suppress component v at 
all nodes in the x<9z-plane and component u at point b. The 
circles in figure 2(b) denote comer nodes and the asterisks 
denote midside nodes. Midside nodes are not present in 
discretizations with four-node tetrahedral elements and eight- 
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z, w 



Tetrahedrons 

Nodes 

1 

1,2, 3, 7 

2 

1,2, 7,5 

3 

2, 7, 5, 6 

4 

1,3, 4,8 

5 

1,7, 8, 5 

6 

1,3, 8,7 



Figure 2.— Finite element models of a cantilever beam. 

(a) Geometric characteristics and loading of beam. 

(b) Finite element discretization, (c) Division of a hexahedral 
ceil into six tetrahedrons. 


node hexahedral elements. The dashed lines in figure 2(b) 
denote hexahedral elements of distorted shapes.Discre- 
tizations with tetrahedral elements were performed by dividing 
each hexahedron into six tetrahedrons, as suggested in 
reference 1. Figure 2(c) shows a typical eight-node hexa- 
hedron divided into six four-node tetrahedrons. The table in 
this figure defines the connectivities of the resulting tetrahe- 
dral elements. A hexahedral volume was divided into six 
10-node tetrahedrons by introducing additional nodes in the 
centers of the edges of the original 4-node tetrahedrons. 

Let us first consider the displacements of the beam. The 
convergence of displacement w at the free end was studied. 
The results are shown in figure 3 for tetrahedral elements, 
and in figure 4 for hexahedral elements. Also shown in fig- 
ures 3 and 4 are values obtained by using corresponding 
isoparametric displacement-based elements. The displacement 
of the free end was normalized with a closed-form solution 
w exact , which was obtained from one-dimensional beam theory 
and includes the effect of shear stresses. Figure 3 shows that 
with 10-node tetrahedral elements, fast convergence is 
achieved, whereas with 4-node tetrahedral elements, conver- 
gence is very slow. Figure 3 also shows that accuracy is not 
improved by increasing the number of independent forces in 
stress interpolation polynomials for four-node tetrahedrons. 



Figure 3. — Convergence of tip displacement of cantilever 
beam studied by using tetrahedral elements. 


These observations agree with those made in the analysis of a 
similar problem in which two-dimensional finite element 
discretizations (ref. 23) using 6-node and 3-node triangles 
correspond to 10-node and 4-node tetrahedrons, respectively. 

The convergence study using hexahedral elements was first 
performed for elements of rectangular shape. The results of 
using the present elements, together with those obtained with 
the standard displacement method, are given in figure 4(a). 
All 20-node elements from the present library provide very 
accurate results with a relatively small number of indepen- 
dent forces. A fast convergence is also achieved with element 
HEX08_18, which employs incomplete second-order poly- 
nomials in the stress interpolation matrix. Elements 
HEX08_48 and HEX08J33, which, respectively, employ 
complete and reduced quadratic polynomials for the stress 
interpolation matrix, produce stiff models that lead to slower 
convergence. A similar behavior was observed in two- 
dimensional analysis of the beam with a four-node quadrilat- 
eral element with bilinear interpolation of the geometry and 
the displacement fields; in three-dimensional analysis this 
corresponds to eight-node hexahedral elements. Note that in 
figure 4(a), the eight-node isoparametric displacement element 
provides a very slow convergence of displacements, whereas 
element HEX08_18 achieves very good accuracy with a small 
number of independent forces. It can, therefore, be concluded 
that the Integrated Force Method significantly outperforms 
the standard displacement method in the analyses of certain 
classes of problems, such as those involving domains of 
regular shapes. 
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Figure 4. — Convergence of tip displacement of cantilever 
beam studied by using hexahedral elements, (a) Regular 
meshes, (b) Distorted meshes. 


The influence of element shapes on the results was also 
studied by using distorted meshes to achieve convergence. 
Distorted elements can be obtained by moving the comer 
nodes a distance of l m = 0.2 L m , as shown in figure 2(b), where 
L m is the dimension of the regular- shaped element. The 
results are shown in figure 4(b). A significant loss of accu- 



Figure 5. — Stress distribution along line z = z g of cantilever 
beam, (a) Normal stress oy. (b) Shear stress r yz . 


racy is suffered by element HEX08_18, whereas the 20-node 
and 8-node elements with quadratic interpolations of the stress 
field are little affected. Figure 4(b) also demonstrates that 
displacement-based isoparametric elements are more sensi- 
tive to distortion than corresponding Integrated Force Method 
elements. 

The stresses were calculated next for a beam of length 
L = 1 2.0 m and cross section d = H = 1 .0 m; its material prop- 
erties were assumed to be E = 21xl0 7 kN/m 2 and v = 0.3, and 
the intensity of the concentrated load P was 60 kN. The re- 
sults obtained with the present elements at normal stress a y 
and shear stress i yz along the line z = -z g = -0.2887 m are 
given in figure 5, along with those from the displacement 
method and the exact solution from the beam theory. Both 
methods performed well when 20-node elements were used. 
Note, however, that element HEX08_18 significantly outper- 
formed the corresponding displacement-based element, which 
exhibited difficulties in stress predictions similar to those ob- 
served for the displacements. 

Example 2: Pure Bending of a Circular Arch 

A circular arch of radius r Q , clamped at 9 = 0° and sub- 
jected to a concentrated moment of intensity Mat 6 - 90° is 
shown in figure 6(a). The arch is assumed to have a uniform 
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Figure 6. — Circular arch subjected to concentrated 
moment M. (a) One-dimensional model, (b) Three- 
dimensional finite element discretization. 


rectangular cross section of dimensions d by H and to be made 
of a homogeneous and isotropic material with parameters 
E and v. This example is presented to verify the use of the 
present developments in the analysis of objects with curved 
contours. A three-dimensional finite element model of the arch 
is shown in figure 6(b). The asterisks and circles in figure 
6(b) represent midside and comer nodes, respectively, as in 
Example 1. The support conditions for the clamped end were 
modeled similarly to those of the cantilever beam. The con- 
tours of the three-dimensional finite element discretization 
were defined by cylindrical surfaces with radii r a and r b , 
as shown in figure 6(b), where r a = r Q - 0.5 H, and r b = 
r 0 + 0.5 H. The arch was analyzed for r 0 = 1 1 .0 m ,d = 1 .0 m, 
H = 2.0 m, E = 21xl0 7 N/m 2 , v = 0.3, and M = 600 kN/m 2 . 

First, the convergence of the horizontal displacement com- 
ponent u of the free end was studied. The results obtained by 



Figure 7.— Convergence of tip displacement u of circular arch. 


using present hexahedral elements and those obtained by 
using a 20-node isoparametric element with the standard 
displacement method are shown in figure 7. The displacements 
were normalized with respect to the analytical solution given 
in reference 30. The present elements performed well, espe- 
cially the 20-node hexahedrons. The results presented here 
for element HEX08_18 were obtained by using a local coor- 
dinate system with an Ox-axis defined by the element 
centroid O and the centroid of one of the element’s sides; an 
Oy-axis normal to the plane defined by the Ox-axis and the 
centroid of the side adjacent to that defining the Ox-axis; and 
an Oz-axis orthogonal to the Oxy-plane. The results are not 
shown for the case with the local axes parallel to the global 
axes, because oscillations were observed in the response and 
convergence was not achieved. Similar behavior was observed 
in the two-dimensional analysis of this problem with element 
QUA04_05, which may be considered to be analogous to 
element HEX08_18. Both of these elements employ incom- 
plete polynomials with the minimum number of independent 
forces, and both perform extremely well for domains of rec- 
tangular shapes. However, because their stress fields are rep- 
resented by incomplete polynomials, these elements may 
become sensitive to the orientation of the coordinate axes. 
Such behavior demonstrates the benefits of using complete 
polynomials in stress field representations: the resulting 
elements are not sensitive to the orientation of local coordi- 
nate systems. 

Next, the stress distributions for the arch were calculated 
by using 20-node elements. The results for normal stresses a r 
and o t along the line r = r g = 10.423 m are compared in 
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Transverse stress, <r t , kN/m 2 Radial stress, a n kN/m 2 



Figure 8. — Stress distribution along line r = r g of circular 

arch, (a) Radial stress <r r (b) Transverse stress <r f . 

figure 8 with the exact values (ref. 30) and with the results 
from using the 20-node isoparametric element with the dis- 
placement method. A good performance of the present ele- 
ments, especially element HEX2CL90, which incorporates 
complete polynomials, can be seen. There is some loss of ac- 
curacy in the results for a r when elements HEX20_57 and 
HEX20__60 are used; this is attributed to using reduced poly- 
nomials in the stress field approximations. The element 
HEX20_90, however, provides better stress predictions than 
the 20-node displacement element. 

£xample3: A Circular Annulus Under Uniform 
Pressure 

A circular annulus with inner radius /?,*, outer radius R 0 , 
and thickness d, subjected to a uniform pressure of intensity 
p along the inner contour, is depicted in figure 9. The annulus 
is assumed to be made of a homogeneous and isotropic mate- 
rial with parameters E and v and assumed to be in the state of 
plane stress. Because of the symmetry of the domain and the 
loading, only a quarter of the annulus was analyzed. A three- 
dimensional finite element model of the annulus is shown in 
figure 9(b). The part of the domain analyzed was discretized 
with two 20-node hexahedral elements in the circumferential 
direction and N elements in the radial direction. The state of 
plane stress in the finite element model was achieved by sup- 



(a) 



Figure 9.— Annular disk under uniform pressure, (a) Two- 
dimensional model, (b) Three-dimensional finite element 
discretization. 


pressing the w displacement component of all nodes lying in 
the plane z = -0.5 d. The symmetry boundary conditions were 
modeled such that u = 0 for nodes in the x = 0 plane, and v = 0 
for nodes in the y = 0 plane. Numerical values for the param- 
eters of the annulus were taken as R t = 5.0 m, R 0 = 10.0 m, 
d- 1.0 m, £ = 21xl0 7 kN/m 2 , v = 0.3, and p = 10 kN/m 2 . A 
convergence study of the radial displacement u of the inner 
and the outer surface was performed first. The results obtained 
with the present elements and those obtained with the 20-node 
isoparametric displacement method element are shown in fig- 
ure 10. A fast convergence of the results can be seen. 

Stress distributions along the radius of the annulus were 
obtained next for N ~ 5 elements in the radial direction. The 
distribution for normal stresses O r and G ti as determined with 
the present elements and with the displacement method, are 
compared in figure 11 with corresponding analytical solutions 
(ref. 30). Again, there is good agreement of the results. Note 
that the locations where stresses were calculated correspond 
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Figure 1 0. — Convergence of annulus radial displacement 
studied by using hexahedrai elements, (a) Inner surface, 
(b) Outer surface. 




(b) Radial distance, r, m 

Figure 1 1 .—Stress distribution along radius of annulus, 
(a) Transverse stress <r t . (b) Fladial stress cr r 



to Gauss integration points for the 2 x 2 x 2 rale. These points 
have been shown to be optimal sampling points in the dis- 
placement and hybrid methods (ref. 9). The stresses were also 
calculated at 3 x 3 x 3 Gauss locations, with no loss of accu- 
racy. Thus, it may be concluded that the Integrated Force 
Method provides a better overall stress response than the dis- 
placement method, which generally predicts stress less accu- 
rately at locations other than the optimal points. 

Example 4: A Simply Supported Rectangular Plate in 
Bending 

A rectangular plate of side lengths 2 a and 2b and thickness 
h is shown in figure 12. The plate is simply supported along 


all four sides and is subjected to a uniformly distributed load 
of intensity p in the direction of the positive Oz-axis. The 
material of the plate is assumed to be homogeneous and iso- 
tropic, with parameters E and v. Because of the symmetry of 
its geometric properties and loading, only a quarter of the plate 
was modeled by using three-dimensional finite elements (see 
fig. 12). The response of the plate was obtained with both 
8-node and 20-node elements. A 4 x 4 mesh was used for 
discretization with the 20-node elements, and a 6 x 6 mesh 
was used for the 8-node elements. Numerical values for geo- 
metric and material parameters of the plate and the loading 
were taken to be a = b = 20.0 m, h = 0.5 m, E = 21xl0 7 
kN/m 2 , V = 0.3, and p- 1.0 kN/m 2 . The stress distribution for 
the plate was calculated first. The stresses <3 X and <5^ along 
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Element type 


□ HEX20JJ0 



(a) Distance, x, m 



Figure 13.— Stress distribution in plate along liney = x. 
(a) Normal stress <r x . (b) Shear stress x^. 


the line y = x and for z = z g = 0.1057 m were calculated by 
using the present elements. They are shown in figure 13 to- 
gether with the stresses obtained by using the standard dis- 
placement method with a 20-node isoparametric element. A 
comparison with the analytical solution given in reference 3 1 
shows good agreement. The displacements of the plate were 
calculated next. Lateral displacements for locations along the 
line x-a were calculated by using both an 8-node and a 20- 
node element (see fig. 14). All 20-node elements, as well as 
the element HEX08_18, provided very accurate displacement 
predictions, but elements HEX08_33 and HEX08_48 pro- 
duced overly rigid models. The element HEX08_18 again 
yielded significantly better results than the corresponding 
eight-node isoparametric displacement element, thereby con- 
firming the conclusions drawn from the results in Example 1. 



2.8x10-3 



Figure 14. — Displacement distribution in plate along line 
x = a. (a) Eight-node elements, (b) Twenty-node elements. 
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Concluding Remarks 

The Integrated Force Method was extended in this paper to 
three-dimensional structural analysis. A general formulation 
was developed to generate stress field interpolations in terms 
of complete polynomials of the required order. Such derived 
stress fields identically satisfied Navier’s equations of equi- 
librium, and the resulting element matrices were invariant with 
respect to the orientation of local coordinate axes and free of 
spurious zero-energy modes. The effect of reducing the num- 
ber of independent forces was also studied. The stress fields 
were derived in terms of reduced polynomials and then, 
expressed in terms of complete and reduced polynomials , were 
used to develop a comprehensive finite element library for 
three-dimensional structural analysis by the Integrated Force 
Method. Both tetrahedral- and hexahedral-shaped elements 
capable of modeling domains with arbitrary configurations 
were developed. To assess the validity and accuracy of the 
elements and to compare the Integrated Force Method with 
the standard displacement method, a number of example prob- 
lems, whose analytical solutions were known, were solved 
with the developments presented herein. The following ob- 
servations can be made on the basis of these numerical 
experiments: 

1. Good accuracy was achieved with all 10-node tetrahe- 
dral and 20-node hexahedral elements. 

2. The elements that employed the complete polynomials 
in stress interpolations exhibited the best overall performance 
and reliability. 

3. Although reducing the polynomials used in stress field 
approximations had no effect on the performance of the cor- 
responding elements in rectangular domains, a certain loss of 
accuracy was observed in the analysis of domains with curved 
boundaries. 


4. The element HEX08_18 performed extremely well in 
analyzing rectangular-shaped objects. Applying the element 
HEX08_18 in the analysis of domains with curved bound- 
aries revealed its sensitivity to the orientation of local coordi- 
nate axes. 

5. Although good results were obtained with a set of local 
coordinate systems that followed the curvature of the domain 
boundaries, unreliable predictions were obtained when the 
local systems were aligned parallel to the global coordinate 
system. Thus, element HEX08_18 should be used with 
caution in the analysis of domains with curved boundaries. 
Such behavior by element HEX08_1 8 also justifies the imple- 
mentation of complete polynomials in stress interpolation 
matrices. 

6. Comparisons of the Integrated Force Method and the 
standard displacement method revealed good performances 
by both methods when 10-node tetrahedrons and 20-node 
hexahedrons were used. 

In most cases, the Integrated Force Method performed bet- 
ter overall in stress calculations. However, the eight-node 
Integrated Force Method elements from the present library 
demonstrated superior behavior in comparison to the corre- 
sponding displacement method element. Thus, for certain 
classes of problems the Integrated Force Method proved to 
be the preferred method of analysis. In general, the Integrated 
Force Method can serve as an alternative to other available 
formulations. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, June 1, 1995 
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Appendix A 

Expressions for Stress Fields 

(a) Complete Polynomials Derived Using the Stress Functions: 

- constant terms: 

4 o) = F} 0) (Ala) 

4 0) = (Alb) 

4 o) = 4 0) (Ale) 

rW = F[ o) (Aid) 

rff = 4 0) (Ale) 

r£> = F 6 (e) (Alf) 

- linear terms: 

4 1} = F^x + F^y + F^z (A2a) 

4 1] = F^x + F^y + F^z (A2b) 

4 1] = F^x + F^y + F^z (A2c) 

rW = F^px — (F^ + Ffp)y + F$z (A2d) 

rff = 4 ) * + F2 ) y-(Fi 1) + F 1 ( 1 1 ^ (A2e) 

= -(F 9 (1) + Fff)x + F$y + F^z (A2f) 

- quadratic terms: 

4 2) = fJV + F 2 (2) y 2 + F 3 (2) z 2 + F 4 (2) xy + F 5 (2) yz + F^zz (A3a) 

4 2 > = F^x 2 + Fj 2) y 2 + F 9 (2) 2 2 + F$xy + F^yz + F$zx (A3b) 

4 2) = F$x 2 + F< 2) y 2 + F x (2) z 2 + F^xy + F$yz + F$zx (A3c) 

ri 2) • = f} 2) ® 2 + F 2 (2) y 2 + F^z 2 - (F< 2) + Fj 2) - F$)xy - (F< 2) + 2F 2 (2) )yz 

- (F{ 2) + 2F 2 (2) )z® (A3d) 

rg> = F 22 } x 2 + F 2 2) y 2 + F 2 (2) z 2 - (F<? + 2F 2 (2) )xy + (F< 2) - fJ 2) - F< 2) )yz 

- (F$ + 2F< 2) )zz (A3e) 

rj 2 > = f| 2 ^® 2 4- Fjg } y 2 4- F$z 2 — (F$ 4- 2Fj^)®y — (F| 2 ^ 4- 2F^)yz 

- (F^-F^ + Fifjyz (A3f) 

- cubic terms: 

4 3) = F 1 (3) x 3 + F 2 (3) y 3 + F 3 (3) ® 3 +Fi 3) ® 2 y + F| 3) ®y 2 

4- F 6 (3) y 2 z 4- F^yz 2 + F^ 3) z 2 x + F^zx 2 - 2 (fJ 3) + F ( $)xyz (A4a) 
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aj 3) = F< 3 ) x 3 + F^y 3 + fJ 3 ) z 3 + F^x 2 y + F$xy 2 

+ Flfy 2 z + F$yz 2 + P$? 9 + F$zx 2 - 2 (F 3 ? + F§)xyz (A4b) 

<rf) = F$x 3 + F$y 3 + F^z 3 + F 2 ( 3 ) x 2 y + F^xy 2 

+ Fify 2 z + F^yz 2 + F 2 ( 3 ) z 2 x + F$zx 2 - 2(F 3 (3) + F$)xyz (A4c) 

4f = F$x 3 + F^y 3 + F$z 3 4- ^(F 2 (3) - 3F& - F< 3 ) )x 2 y 

+ |(F# - 3F<? - F 4 (3) )xy 2 + F$y 2 z- (if } + ZF$)yz 2 - (F$ + 3 F$)z 2 x 
+ Fg>zx 2 + (3F 2 (3) - F 9 (3) - F$)xyz (A4d) 

rg> = Fg } x 3 + F^y 3 + F 3 {3) * 3 - (Fg> + 3Fg>)x 2 y + F^xy 2 

4- i(F 4 (3) - 3 F<? - F^)y 2 z - |(F 9 (3) - 3F 2 (3) - F$)yz 2 + F 3 (3) * 2 x 

- (F< 3) + 3 F$)zx 2 4 - (3F< 3) - F<? - F$)xyz (A4e) 

rj 3 > = F$x 3 + F 3 ( 3 ) y 3 + F$z 3 + Fjg>x 2 y - (F 2 (3) 4- 3F^)xy 2 

- (F 5 (3) 4- 3F^)y 2 z - Fg>yz 2 + |(^ 3) - 3 F 2 ( ? - fJ 3 ) )z 2 x + 

+ \(fS } ~ 3F} 3) - F$)zx 2 + (3F} 3) - F 4 (3) - F$)xyz (A4f) 


The superscripts 0, 1, 2, and 3 in Eqs.(Al) - (A4) denote the stress components and the independent 
forces that correspond to constant, linear, quadratic and cubic terms in the stress interpolations, 
respectively. 

- complete cubic polynomial: 

cr s = Fi 4- F 7 X + Fsy + F 9 z 4- F 22 x 2 4- Fwy 2 4- F 24 z 2 4- F 25 xy + F 26 yz 4- F 27 zx 

+ F49X 3 + F 50 y 3 + F51Z 3 4- F 52 x 2 y + F 53 xy 2 + F 54 y 2 z + F S5 yz 2 + F 56 z 2 x 4- F 57 zx 2 

- 2(F 7 q + F$o)xyz (A5a) 

<Ty — F 2 + F10X 4- F\\y 4- Fl 2 z 4- F 23 x 2 4- F 2 gy 2 + F^qz 2 + Fzixy + F 32 yz 4- F 33 zx 

4* F 58 x 3 4- Fsgy 3 4- Fgoz 3 4- F 61 x 2 y + Fe 2 xy 2 4- F^y 2 z 4- F 64 yz 2 4- F 6 sz 2 x 4- Fqqzx 2 

- 2(Fso 4- F S s)xyz (A5b) 

<r* = F 3 4- Fi 3 x 4- Fi 4 y 4- F\ 9 z 4- F 34 X 2 4- F 3 gy 2 4- Fz&z 2 4- F 3 rxy 4- F&yz + F&zx 

4- F 6 7 X 3 + Fs$y 3 + F 69 z 3 4- F 70 x 2 y + F 7 l xy 2 + F 72 y 2 z 4- F 73 yz 2 + F 74 z 2 x 4- F 75 zx 2 

- 2(F§4 4- Fs 9 )xyz (A5c) 

T xy = F 4 + F 17 x - (F 7 4- F 2 i)y + F 16 z 4- F 40 x 2 4- F 41 y 2 4- F 42 z 2 

- (F 22 4- F 2 g — F 3 6)xy — (F 27 4“ 2 F 49 )yz — (F 32 4- 2F 4 s)zx 

4- F 7 gx 3 4- F 77 y 3 4- F 73 z 3 4- 2 (^ 74 — ^F 4 9 — F 6 2 )x 2 y 4- -(F 73 — 3 F 59 — Fs 2 )xy 2 
4- F 79 y 2 z - (F 56 4- 3Fss)yz 2 - (F^ 4- 3F 83 )z 2 x 4- F^zx 2 

4- (3F 69 — F5 7 — F&i)xyz (A5d) 
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T yz = Fs + F\sx + F\sy — {Fii + Fn)z + F&x 2 + F 44 y 2 + F45Z 2 

— (F39 + 2F4G,)xy + (F22 — F29 — F^)yz — (F3x + 2 F 4 o)zx 

+ F 81 x* 4 - F 8 2y 3 + F&z 3 — (F75 + ZFse)x 2 y + F&^xy 2 + -(F52 — 3F59 — Frz)y 2 z 

— 2 ^ 57 ~ ^69 — Fsa)yz 2 + Fg$z 2 x — (F$i + 3 F 76 )zx 2 

+ ( 3 F 49 - F 62 - F 74 )xyz (A 5 e) 

t zx = Fq — (F15 + Fig)x + F 20 y 4 - F 2 xz + F 4 6X 2 4- F 47 y 2 4- F 4 s.z 2 

— {F33 + 2 F 44 )xy — (F25 + 2 F 4 i)yz - (F22 - F 2 9 + Fz&)zx 

4 - F%$x 3 + F 87 y 3 + Fgsz 3 + F 88 x 2 y — (F72 + 3 Fs 2 )xy 2 — (F53 4 - ZF 77 )y 2 z 4- Fgoyz 2 
+ 2^ 63 — ^ 69 — -^ 57 )^ 2 * + 2(^62 — 3 F 4 g - F 7 4 )zx 2 
+ (3F59 — F52 — F 78 )xyz 

(b) Robinson’s Stress Field [5]; 

o x = Fi + F 2 y 4 - F 3 z + F 4 yz 
o y = F 5 + F 6 z 4 F 7 x + F 8 zx 
o z = Fg + F 10 X + Fny 4 Fi 2 xy 
Tiv = 1*13 + F 14 Z 

Tyz = F15 + F\qX 

T zx = -F17 4 - Fisy 

(c) Spilker’s Stress Field [9]: 

Ox — Fj — (F13 4 - F 24 4 - 2 F 36 + 2F5i)x + F 2 y 4 - F$z 

— (F 7 + Fg + F18 4 - F 2 g)x 2 + 2 F 7 y 2 + 2 F 8 z 2 
4 - (F 4 — F15 — F 2 6 — 2 F 4 x)xy + Fsyz + (Fg — F17 — F 2 s — F39)zx 
+ g ( 2 F 22 4 - 2F33 + F54 + Fse)x 3 

— j (F10 4 - Fa + F 2 x — F 22 4 - F3 2 — F33 + F5 2 4 F53 — F5 4 4 F55 — F56 + F57) (y 3 4 z 3 ) 

4 F 52 x 2 y - F&xy 2 + Fx 0 y 2 z + Fnyz 2 - F 22 z 2 x 4 F 53 zx 2 + (F 9 - F 20 - F 3 x)xyz (A 7 a) 
o y = Fx 2 4 F13X — (F 2 + F 2 5 4 2F35 4 2 F 4 6 )y + Fi 4 z 
4 2 Fl 8 X 2 — (F7 4 Fi8 + F19 4 F3o)y 2 4 2 Fiqz 2 

— (F 4 — F15 4 F 2 6 4 2F 40 )xy — (F5 — Fi6 + F 27 + 2 F 42 )yz + F\ 7 zx 
+ g( 2 Fn 4 2F3 2 + F5 2 4 F 87 )y 3 

— g ( 1*10 - 1*11 + 1*21 + F22 ~ 1*32 + 1*33 — F52 4 F53 4 F5 4 4 F55 4 F56 — F57) (z 3 + X 3 ) 

— F3 2 x 2 y 4 F$4xy 2 4 Fssy 2 z — Fnyz 2 4 F 22 z 2 x 4 F 2 xzx 2 


(A 5 f) 


(A6a) 

(A6b) 

(A6c) 

(A6d) 

(A6e) 

(A6f) 
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— (Fg — F 2 o + Fz\)xyz (A7b) 

&z = i *23 + F 24 x + F 25 y - (F3 + F14 + 2F45 + 2 Fso)z 

+ 2 F 2 gX 2 + 2F30J/ 2 — (Fg + Fig + F29 + F$o)z 2 
+ F 26 xy — (Fg + Fi6 — F 2 7 + 2 F 4 s)yz — (Fg + Fn — F 2 8 + 2F 4 j)zx 
+ g (2Fio + 2F 2 1 + F53 + Fgg)z Z 

+ g(Fio — Fn + F 2 x — F 22 — F3 2 — F33 — Fg 2 + F53 — F54 + F55 — F56 — F57) (® 3 + y 3 ) 

+ F 32 x 2 y + Fzzxy 2 - F w y 2 z + F$ 7 yz 2 + F 5 Q z 2 x - F 2 \zx 2 

— (Fg + F 20 - Fzi)xyz (A 7 c) 

T X y — F34 4- (F 2 + 2 F 3 g)x + (F13 + 2 F 3 e)y 4- (F37 + F3s)z 

+ jW - jfis)* 2 - 5(5*1 - W + Wo + F«)z ! 

+ 2(F7 + Fis)xy — (-Fe — Fn — 2 F 3 g)yz -f (F5 — -F16 4 - 2 F 42 )zx 

— ^(Fio + F 2 i + F53 + Fgg)z 3 — (F 22 + F33 + F54)x 2 y — (Fn + F3 2 + Fg 2 )xy 2 

— \(\f* - F 20 )y 2 z + F 22 yz 2 + F n z 2 x + i(F„ - \f 2Q )x 2 z 

+ 2 (Fio + F 2 i)xyz (A7d) 

Tyz = F43 + (F38 + F44)* + (F14 + 2 F 4 g)y + (F 2 5 + 2 F io)Z 
+ (F42 + F 4 s)x 2 + -(Fi6 — ^F 27 )y 2 - -(-Fi6 ~ 

+ (Fn — -F 2 s + 2 F 47 )xy + 2 (Fig + Fzo)yz - (~F \ 5 - F 2 6 — 2 F 4 o)zx 

— g (F 22 + F33 + F54 + F56)* 3 + F 2 \x 2 y + ~(F 2 o — g F3i)a:y 2 + (F10 + F 2 i — Fss)yz 2 

— (Fn + F3 2 + Fg 7 )y 2 z — -(-F 2 o — F3i)z 2 x + F 32 zx 2 

+ 2 (F 22 + Fzz)xyz (A 7 e) 

r z x — F49 + (F3 4- 2Fgo)x 4- (F37 + F 44 )y 4 - (F 2 4 + 2Fg\)z 
+ |(F 6 - \f 2S )x 2 + (F 39 + F 47 )y 2 - |(|f 6 - F 28 )z 2 

+ (F5 - -F 2 7 + 2 F 48 )xy - (~F 4 — F 2 6 — 2 F 4i )yz 4 - 2 (Fs 4- F 2 g)zx 

— + ^32 4- Fg 2 4- F57)y 3 4- -(Fg — -Fzi)x 2 y + Fioxy 2 

+ Fzzy 2 z — -(-F 9 — Fz\)yz 2 — (F 22 4- F33 4 - Fg&)x 2 z — (F10 4- F 2 i + F 83 )xz 2 

4- 2(Fn 4- F 32 )xyz (A7f) 

(e) Quadratic Terms for the Stress Field in the Element HEX08.33: 

4 2r > = F 1 (2r) y 2 + F 2 (2r) z 2 + F 3 (2r) yz 
4 2r ) = Fi 2 r) z 2 + Ft ) x 2 + Ft ) zx 
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(A8a) 

(A8b) 



a (2r) = F j?r) x 2 +F Vr ) y 2 + F {2r) xy (A8c) 

= F$ r) z* (A8d) 

4l r) = 4 2r) * 2 (A8e) 

4l r) = HPy 2 (A8f) 

(f) Quadratic Terms for the Stress Field in the Element THD10.39: 

a (2r) = p(2r) x 2 + p (2r) y 2 + p {2r) j + F }?%y + F & r ) yz + pi*) zx (A9a) 

c$ r) = F? r) x 2 + F 8 (2r) t/ 2 + Fg 2r) z 2 + F{o r) xy + Ftf r) yz + F% r) zx (A9b) 

4 2r) = F^ ) x 2 + F^ r) y 2 + F^ ) z 2 + F^ ) xy + F! 2r) yz + F^ ) zx (A9c) 

rg r) = (-if r) - F 8 (2r) + F^ r) )xy - F 6 (2r) yz - F^zx (A9d) 

rg r) = -Fj.l r) xy + (ff r) - F^ r) - F^)yz - F^zx (A9e) 

r< 2r ) = -F$ r) xy - F {2r) yz - (F} 2r) - F 8 (2r) + F< 2r) )z* (A9f) 

(f ) Quadratic Terms for the Stress Field in the Element TET10.36: 

4 2r) = Ff 2r) xy + F? r) yz + F^zx + F x (2r) x 2 + F^ r) {y 2 + z 2 ) (AlOa) 

4 2r) = Fj 2r) xy + Fi 2r) yz + Ff r) zx + F{ 2r) y 2 + F< 2r) (z 2 + x 2 ) (AlOb) 

4 2r) = Fi 2r) xy + F^yz + F^zx + F^z 2 + F< 2r) (x 2 + y 2 ) (AlOc) 

4 2r) = (-F§ r) -F^P + F$ ) )xy-Ft ) yz-Fi 2r \zx (AlOd) 

rg r) = -F 9 (2r) *y + (Fg r) - F<? r) - F$ r) )yz - f[ 2t) zx (AlOe) 

4 2r) = -Fi 2r) xy-F} 2r) yz+(Ffi r) - F{ 2r) + F{ 2r) )zx (AlOf) 

(h) Cubic Terms for the Stress Field in the Element HEX20.60: 

4 3r) = F} 3r) yz 2 + pf r) y 2 z + F 8 (3r) s 3 + F$ 3r) xyz + ZF$ r) x 2 z (Alla) 

4 3r) = Fi 3r) z 2 x + Fj? r) zx 2 + F§ r) y 3 + F^xyz + ZF^xy 2 (Allb) 

4 3r) = F^xy 2 -f F 6 (3r) x 2 y + F^z z + F^xyz + ZF^yz 2 (Allc) 

4 2r) = -3F s x 2 y - ^F 7 y 2 z (Alld) 

= -ZF 10 y 2 z - i F g xz 2 (Alle) 

r 2 (3r) = -ZF u z 2 x - ^F n xy 2 (Allf) 
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Appendix B 
Symbols 


{»,} 

element equilibrium matrix; J\z) T [Y]dV 

Ns 

number of prescribed displacements 

[B p ] 

portion of system equilibrium matrix corre- 
sponding to nodes where external loads are 

Nt 

total number of displacement degrees of 
freedom 


prescribed 

n 

number of system equilibrium equations; 

[Bj] 

portion of system equilibrium matrix corre- 


— N s 


sponding to nodes with prescribed displacement 
boundary conditions 

Ox, Oy, Oz 

global coordinate axes 

[C] 

compatibility matrix 

{P} 

vector of system equivalent nodal loads 

c*. 

arbitrary coefficients 

{P e } 

vector of element equivalent nodal loads 

[D] 

compliance matrix of material 

[P*] 

ii 

Si 

d 

cross sectional dimension 


L 1 J J 

E 

modulus of elasticity 

P 

order of polynomial 

{F} 

system vector of independent forces 

Q 

number of independent generalized forces 

{F e } 

vector of element independent generalized 

{R} 

vector of support reactions; {B 5 ]{F} 


forces 

Ri>Ro 

inner, outer radius of circular annulus 

F; 

generalized force coefficients 

r 

number of compatibility conditions 

[G] 

system flexibility matrix 

r,9 

polar coordinates 

{G e } 

element flexibility matrix; f (Y) T [D][Y]dV 

jy 

r a 

radius to inside of arch 

H 

cross sectional dimension 

r b 

radius to outside of arch 

[J] 

n x m deformation matrix; represents top n rows 

oftsr 1 

r 8 

radius of Gauss point 



r 0 

radius of circular arch 

[L] 

matrix of differential operators that defines 
strain-displacement relationship 

length 

[S] 

[s KSJ 

L 


lm 

distance comer node is moved 

{U} 

vector of unknown nodal displacements; 
[J][G]{F} 

M 

intensity of moment 

{U e } 

vector of displacement at element nodes result- 

m 

number of independent forces 


ing from finite element discretization 

[N] 

matrix of displacement interpolation functions 

{«} 

displacement vector 
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w 

displacement at free end of beam 

V 

Poisson’s ratio 

^exact 

exact value of displacement 

W T 

stress vector of {a x o y O z x^ x yz x^} 

[Y] 

stress interpolation matrix 

Gr&t 

radial and tangential stress components 

[Z] 

MM 

Gx&y>®z 

normal stress components 

z g 

^-coordinate of Gauss point 

x xy>Tyz> x zx 

shear stress components 


vector of element deformations corresponding 
to forces {F e } 


stress functions 




linear functions of constants 
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